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1.  Introduction 

The  approach  to  seismic  source  discrimination  through  the  joint  use  of  seismic,  direct 
acoustic  and  indirect  atmospheric  wave  sensing  using  electromagnetic  measurements  (passive 
or  active)  is  illustrated  in  Table  1.  Here  seismic  events  of  different  types,  made  comparable 
by  normalizing  to  the  same  low  frequency  compressional  wave  amplitude,  are  compared  in 
terms  of  qualitative  estimates  of  the  amplitude  of  the  different  diagnostic  signals  relative  to 
noise  background.  These  relative  signal  amplitudes  are  indicated  in  the  columns  under  each 
source  type  as  being  "large",  "moderate"  or  "small"  relative  to  the  typical  noise  level.  These 
estimates  are  based  on  seismic  observations  ( e.g .,  Evemden  et.  al„  1986)  and  ionospheric 
acoustic  wave  measurements  using  EM  sensing  (e.g.,  Blanc,  1982),  together  with  rough  esti¬ 
mates  for  excitation  of  secondary  ionospheric  EM  emissions  and  the  direct  near  surface  acous¬ 
tic  or  gravity  wave.  The  expected  "signatures"  of  the  sources  in  terms  of  the  rough  sizes  of 
the  various  signals  are  also  indicated,  along  with  the  basis  for  discriminating  between  the 
different  source  types. 

In  order  to  assess  quantitative  signal  levels,  we  have  developed  and  applied  atmospheric 
modeling  methods  to  provide  a  basis  for  discrimination  of  seismic  events  using  a  combination 
of  electromagnetic  and  seismic  sensing  methods  to  identify  small  chemical  and  nuclear  tests, 
as  well  as  earthquakes.  The  primary  objective  of  the  research  has  been  to  predict  low  fre¬ 
quency  gravity  waves  in  the  atmosphere,  produced  by  surface  and  buried  explosions,  that  pro¬ 
pagate  to  high  altitudes  and  produce  large  amplitude  waves  in  the  ionosphere.  Since  these 
waves,  which  increase  in  amplitude  with  altitude  because  of  the  decreasing  density  of  the 
atmosphere  with  height,  will  produce  relatively  large  fluctuations  in  the  electron  densities  in 
the  ionosphere,  the  disturbance  can  be  sensed  by  standard  electromagnetic  sounding  tech- 
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niques.  Hence  sensitive  monitoring  of  explosion  produced  atmospheric  disturbances  can  be 
accomplished  by  EM  sounding  methods.  Our  objective  is  therefore  to  provide  predictions  of 
the  ionospheric  disturbances  to  be  expected  from  different  kinds  of  shallow  seismic  sources  so 
that  the  characteristic  atmospheric  wave  signatures  of  these  sources  can  be  used,  along  with 
seismic  methods,  to  help  identify  them.  We  also  use  the  current  nonlinear  atmospheric 
modeling  capability  to  systematically  study  coupling  between  near  surface  atmospheric  tur¬ 
bulence  and  seismic  noise,  in  a  variety  of  earth  models,  in  order  to  more  fully  understand  and 
predict  high  frequency  seismic  noise  variability  in  different  geologic  environments. 

Another  principal  objective  of  this  research  has  been  to  provide  predictions  of  complete 
local  and  regional  seismic  wave  fields  to  be  expected  from  explosions  (single  and  multiple) 
and  shallow  earthquakes  in  complex,  laterally  varying  two  and  three  dimensional  anelastic 
earth  models.  In  this  work,  particular  importance  was  placed  on  achieving  predictive  capabili¬ 
ties  over  the  seismic  frequency  band  from  .02 Hz  to  5Hz,  where  representation  to  relatively 
high  frequencies,  at  or  above  5Hz,  was  considered  to  be  of  major  importance  in  view  of  the 
low  magnitude  events  that  are  of  primary  interest  for  regional  seismic  discrimination.  On  the 
other  hand,  the  lower  frequency  end  of  the  range  is  important  in  the  prediction  of  atmospheric 
coupling,  particularly  for  estimating  the  excitation  of  low  frequency  gravity  waves. 

Achieving  results  that  are  close  to  the  observed  complexity  of  seismic  wave  fields, 
through  incorporation  of  vertical  and  lateral  randomness,  rough  topography  and  abrupt  lateral 
changes  in  crust-upper  mantle  seismic  velocity  structure  was  considered  important  for  seismic 
discrimination,  with  some  of  this  medium  variability  being  important  in  characterizing  atmos¬ 
pheric  wave  excitation  as  well.  To  obtain  realistic  estimates  major  effort  involved  develop¬ 
ment  and  application  of  different  2  and  3-D  numerical  modeling  methods  to  quantify  the 
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Discrimination  by  Seismic  Signal(s) 


mechanisms  of  generation  of  discriminatory  signals  and  to  quantify  their  variability  as  func¬ 
tions  of  changes  in  3-D  earth  structure  and  source  type. 

2.  Modeling  Seismic  Signal  Propagation  in  the  Near  and  Regional  Distance  Ranges 

2.1.  Numerical  Methods  for  Seismic  Wavefield  Modeling 

Most  modeling  of  seismic  wavefields  in  laterally  heterogeneous  media  has  used  spatially 
discrete  numerical  methods  such  as  the  finite  difference  or  finite  element  technique.  The 
main  reason  for  using  a  pseudospectral  method  instead  of  these  conventional  techniques 
is  the  increased  accuracy  of  a  pseudospectral  approximation,  which  allows  for  larger  scale 
and/or  higher  frequency  simulations.  The  pseudospectral  method  actually  belongs  to  a 
larger  class  of  discretization  techniques  known  as  spectral  methods.  The  significance  of  the 
nomenclature  is  explained  below.  Comparisons  have  been  made  between  finite  difference  and 
Fourier  pseudospectral  methods  in  terms  of  runtimes,  memory  requirements  and  accuracy  of 
solutions  [6]  [3]  [14].  The  general  conclusion  is  that  for  two-dimensional  simulations  of  the 
same  wavefield,  a  fourth-order  finite  difference  method  currently  runs  about  twice  as  fast  as 
a  Fourier  method,  but  requires  about  twice  the  memory.  Of  course  a  specific  comparison 
between  meti  ods  requires  the  consideration  of  how  they  incorporate  the  following  factors, 
among  others:  sources,  absorbing  boundaries,  anelastic  attenuation,  material  structure  sam¬ 
pling,  a  free  surface  condition,  and  machine  vectorization  and  concurrency.  Although  no  firm 
conclusions  have  been  drawn  for  three-dimensional  problems,  storage  requirements  alone  fa¬ 
vor  the  use  of  spectral  methods  for  large  calculations.  This  has  led  us  to  investigate  their 
usefulness  for  both  2-D  and  3-D  modeling. 

The  spectral  methods  we  are  investigating  solve  the  elastodynamic  equations  of  motion 
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by  approximating  the  spatial  dependencies  as  truncated  series  of  orthogonal  functions  and 
by  integrating  the  expansion  coefficients  in  time  as  in  a  finite  difference  method.  We  derive 
the  method  from  a  variational  formulation  of  momentum  conservation,  and  we  include  a  free 
surface  boundary  condition  by  making  the  spatial  domain  nonperiodic.  Although  solving 
nonperiodic  problems  is  usually  computationally  much  more  costly  than  periodic  problems 
that  use  harmonic  functions,  our  method  requires  little  more  computation  than  the  Fourier 
method,  and  it  accurately  simulates  surface  waves. 

2.2.  Variational  Formulation  of  Spectral  Methods 

In  the  equations  that  follow,  Greek  subscripts  denote  spatial  coordinate  directions,  and 
tip  is  a  unit  vector  in  the  direction.  Boldface  subscripts  and  superscripts  on  summations 
represent  three-dimensional  sets  of  integers,  e.g.  k  =  k3),  so  that  represents  a 

triple  sum.  Summation  over  repeated  indices  is  assumed,  and  the  symbol  *  represents 

For  a  spectral  solution  to  the  elastodynamic  equations  of  motion,  we  expand  each  com¬ 
ponent  of  the  displacement  field  in  a  truncated  series  of  infinitely  differentiable,  orthogonal 

basis  functions  6(k,x)  over  the  volume  Vx  =  FI/Ui  of  the  spatial  domain: 

K /a 

ua(x,t)  =  £  KM)Pa(M)  (1) 

k=— K/2 

Let  the  symbols  (  ),  {  },  and  [  ]  denote  a  row  vector,  column  vector,  and  square  matrix, 
respectively,  so  that  equation  (1)  may  be  written  as 

ua{x,t)  =  <6(x)){<7tt(i)}  (2) 

From  orthogonality,  the  wavenumber  coefficients  are 

{CUt)>  =  L  <6-(x))«„(x, ()<*>*  (3) 
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With  a  strain  field 


«.«(*.<)  =  |  (< £;H*))m‘))  +  ( sf;Kx))(W)))  (4) 

we  consider  a  general  constitutive  relation  and  express  the  stress  tensor  as 

=  Ea0^{x)ta0(x,t)  (5) 

In  order  to  obtain  a  governing  equation  for  the  coefficients  Ua(t),  we  substitute  the  expansions 
of  equations  (2),  (4),  and  (5)  into  the  variational  statement  of  momentum  conservation: 

Jv  [  P(*)  Hju(x,  t)  ■  £u*(x)  +  £(x,  t)  ■  6c*(x)  -  f (x,  t)  -  6u’(x)  J  <fx 

—  i  t(x,  t)  ■  Sum(x)  dS  =  0  (6) 

Js, 

where  ^u(x)  is  a  virtual  displacement,  6e(x)  is  a  virtual  strain,  and  a  *  indicates  a  complex 
conjugate.  /0(x,  t)  is  a  body  force  density  and  the  surface  integral  is  taken  over  that  part  of 


the  surface  on  which  the  tractions  ta(x,  t)  are  applied.  We  express  the  result  as 

{«&.}'(  (*]{  &&(<)>  +  -  {&(<)}  )  =  o  (7) 

where 

(M)  =  /  «x))V(x  M*))#*  (8) 

[*..]  =  Jvj  ^-6(x))<Pi  (9) 

{£.(<)}  =  /  «x))'/,,(x,  t)<Pz  +  /  (4(x))'(.(x.  t)iS  (10) 

JVx  JSx 

and  a  f  represents  conjugate  transpose.  The  null  vector  is  the  only  vector  orthogonal  to  all 
virtual  (  unconstrained  )  displacements,  so  we  have  as  a  governing  equation 

(M]{  £&,(<)}  +  [JU{£>,,(<)}  -  (£,(()}  =0  (11) 
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Before  choosing  a  particular  basis  set  (6(x)},  let  us  change  the  formulation  slightly  in 
anticipation  of  the  numerical  solution.  It  is  convenient  to  keep  the  stress  tensor  explicit  in 
order  to  apply  anelastic  attenuation  in  the  time  domain  as  per  Emmerich  and  Korn  [4)  or 
Witte  and  Richards  [16].  To  that  end,  we  express  the  restoring  force 

{«.(<)}  =  [*.<]{(><«}  (12) 

in  terms  of  the  stress  tensor  by  expanding  the  stress  tensor  with  the  basis  functions: 

^(x,t)  =  (b(x))ftM‘))  (13) 

Expressions  for  the  coefficients  {Ta/j(t)}  in  terms  of  displacement  and  moduli  coefficients  are 
obtained  by  comparing  equation  (13)  to  equation  (5)  once  a  particular  constitutive  relation 
is  chosen.  In  practice,  we  solve  two  first-order  partial  differential  equations  for  stress  and 
velocity  and  write  momentum  conservation  as 

[MH  |v.(l)}  +  {£*]{£*<<)}  -  {£,(<))  =  0  (14) 

fl  *  A 

'  where  {  ffiVa(t)}  is  a  vector  of  velocity  coefficients  and  [Dp]  is  the  divergence  matrix: 

\be\  =  (15) 

Given  a  constitutive  relation  and  an  appropriate  set  of  basis  functions,  we  can  obtain 
explicit  equations  for  the  expansion  coefficients  Va(t)  and  Tap(t).  In  what  follows,  consider 
as  a  special  case  isotropic  material  whose  modulus  is  given  by 

Ea06y(x)  =  A(x)  6a06^  +  p(x)  (6aS6fo  +  S^Ssp)  (16) 
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In  general,  however,  no  restrictions  are  placed  on  the  symmetry  of  the  medium.  The  density 

and  Lame  coefficients  sire  expanded  in  the  chosen  basis,  e.g. 

K/2 

W«)  =  £  E„w,(k)Hx,k)  (17) 

lc=-K/J 

and  the  integrals  in  equations  (8),  (10),  and  (15)  are  solved  by  invoking  orthogonality. 

2.3.  Boundary  Conditions 

Because  equation  (14)  for  momentum  conservation  was  obtained  from  a  variational  princi¬ 
ple,  the  natural  boundary  conditions  on  the  surface  of  the  domain  are  automatically  satisfied. 
In  the  absence  of  the  surface  integral  in  equation  (10),  a  traction-free  boundary  condition 
is  implicit  in  the  formulation,  but  the  basis  functions  themselves  must  allow  the  free  sur¬ 
face  condition.  In  the  Fourier  spectral  method,  the  basis  functions  are  sines  and  cosines 
over  the  interval  2x  in  each  coordinate  direction,  and  therefore  the  boundary  conditions  are 
periodic.  On  the  other  hand,  a  traction-free  boundary  corresponds  to  a  nonperiodic  prob¬ 
lem.  For  nonperiodic  problems,  the  basis  functions  must  be  nonperiodic,  as  are  Chebychev 
or  Legendre  polynomials.  Chebychev  polynomials  have  been  used  successfully  to  solve  the 
free  surface  condition  in  elasticity  [11],  but  the  method  is  not  derived  from  a  variational 
principal.  Using  a  variational  formulation  ensures  that  the  stiffness  matrix  of  equation  (9) 
is  Hermitian  positive  semi-definite  and  hence  the  eigenvalues  of  this  differential  operator  are 
real  and  non-negative.  It  is  well  known  that  spectral  methods  for  non-periodic  domains  lead 
to  non-Hermitian  matrices  [1].  If  one  were  to  use  a  variational  formulation  with  a  Chebychev 
basis  set,  the  integrals  in  equations  (8),  (10),  and  (15)  could  not  be  evaluated  by  invoking 
orthogonality,  and  such  a  formulation  would  be  computationally  prohibitive.  Because  the 
differential  operator  of  a  Chebychev  spectral  method  possesses  complex  eigenvalues,  a  sim- 
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pie  leapfrog  time  integration  scheme  is  not  stable,  and  a  fourth-order  Runge  Kutta  scheme 
is  usually  employed.  The  boundary  conditions  are  made  stable  by  adding  additional  con¬ 
straints  based  on  characteristic  variables  [1],  The  computational  effort  using  this  scheme  for 
a  problem  with  a  given  spatial  domain  size  is  four  times  that  of  a  simple  leapfrog  method  for 
the  same  problem,  and  it  requires  four  levels  of  storage.  Kosloff,  et.  al.  [10]  have  recently 
developed  a  high-order  time  integration  method  referred  to  as  the  Rapid  Expansion  Method, 
which  may  improve  the  efficiency  of  the  Chebychev  method  while  eliminating  numerical  dis¬ 
persion  entirely.  However,  they  note  that  time  histories  must  be  kept  in  storage,  which  seems 
to  be  a  prohibitive  requirement. 

2.4.  The  Fourier  Sawtooth  Method 

Our  approach  to  incorporating  a  free  surface  condition  into  a  spectral  method  has  been  to 
use  a  basis  set  comprised  of  harmonic  functions  (  the  Fourier  method  ),  plus  additional  terms 
in  the  vertical  direction  that  will  decouple  the  ends  of  the  otherwise  periodic  spatial  domain. 
We  use  a  basis  set  that  maintains  orthogonality  so  that  the  scheme  is  computationaly  efficient 
and  can  be  solved  with  a  simple  low-order  time  integration  method.  Polynomials  have  been 
used  with  Fourier  series  to  improve  the  convergence  for  nonperiodic  problems  [8],  but  they 
have  not  been  included  in  a  variational  formulation.  We  consider  here  a  Fourier  set  plus  a 
linear  term.  The  linear  term’s  spatial  dependence  is  that  of  a  sawtooth  minus  the  Fourier 
representation  of  the  sawtooth.  With  this  form  the  sawtooth  function  is  orthogonad  to  all 
Fourier  terms.  Its  spatial  dependence  is  illustrated  in  Figure  (1).  Applying  this  mixed 
basis  set  to  the  vauriational  formulation  of  the  previous  section,  we  derive  expressions  in  two 
dimensions,  but  the  derivation  is  easily  extended  to  3-D. 
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Let  the  field  variables  ’  expansion  be 


K/2  *,/  2 

va(x,t)  =  £  6(k,x)Va(M)  +  £  *(*!,/,*)  V.(*if/,0 

k=-K/2  fc,=-ff,/2 

where 

6(k,x)  =  ei2"kfi*p/xe 


(18) 


(19) 


and 

x)  =  _(A  +  i  j;  X  (20) 

*7*0 

with  normalization 

Wii,/,x)|  =  MiBi  ;  B.  =  ^  -  £  (21) 

ka?to  k 

The  index  I  indicates  terms  associated  with  the  x%  dependence  of  the  sawtooth  minus  the 
Fourier  representation  of  the  sawtooth.  With  this  form  b{kx,  7,x)  is  orthogonal  to  b{ k,x)  for 
all  k.  From  now  on  we  neglect  the  limits  on  summations,  so  that  the  limits  are  implied  to 
be  those  of  equation  (18)  unless  otherwise  indicated.  By  expanding  the  density  and  Lame 
coefficients  in  the  basis  functions  of  equations  (19)  and  (20),  and  by  defining 


**„/)  =  **,,0)  +  *  E  **..*)[  wp  -  rii  £  Ejsfes  ]  (M) 

and  similarly  for  A(&i,/)  and  /*(&i,  /),  equation  (14)  in  the  absence  of  body  forces  becomes 

S«k-D^.(i,0  =  ^r„(k.o+  (23) 

E«k-l)^(1.0  =  ^f1J(k,i)+^fB(k,()  (24) 

=  ^f,1(k„/,0-B^E2'n(l.,l1,<)  (25) 

h  h 

EHh-l„/)$Vi(I„I,t)  =  £t„(k„l,,t)  (26) 

<1  h 


10 


The  stress  tensor  coefficients  are  obtained  from  the  constitutive  relation: 


|r„(k,l)  =  E(  (  A(k  -  I)  +  2/i(k  -  1)  ) 

+  Mk-l)(^Va(I,0+^Vi(J„  /,<))]  (27) 

|r»(k,i)  =  E[Wk-i)^-v2(i,() 

1 

+  Mk  - 1)  (^v,(l,<)  +  £&(/,,  f,t)  )  ]  (28) 

ffWM)  =  E[i0'-l)¥l('1(l,<) 

1 

+  (A(k— l)  +  2/i(k-l))(^»H(l,()+  £v2(l> ,/,!))]  (29) 


(*„/,()  = 

SfTW*.,/,<)  = 
&fn(k„I,t)  = 


E  ( Mk,  -/.,/)+  2A(*. -/„/))^  V,(/„ /, f) 

Jl 

EA(*.-<1,/)T7lfc(/i,/,() 

*1 

E  *<*>  -/>,/) 

1 1 


(30) 

(31) 

(32) 


Notice  that  products  in  the  spatial  domain  have  become  convolutions  in  the  wavenumber 
domain.  Their  general  form  is 


*00  =  ££(k-k)*(  k)  (33) 

]c 

and  they  typically  are  computed  with  a  Fast  Fourier  Transform  (  FFT  ).  The  computational 
efficiency  of  the  FFT  precludes  the  use  of  any  other  (  current )  methods  for  performing  the 
convolutions.  The  wavenumber  coefficients  for  the  velocity  and  stress  fields  are  obtained  by 
numerically  integrating  equations  (23)  -  (32)  in  time. 

25.  The  Collocation  or  Pseudospectral  Method 

If  we  solve  the  governing  equations  in  the  spatial  domain  instead  of  in  the  wavenumber 
domain,  both  domains  are  discretized.  The  two  domains  are  related  by  a  discrete  Fourier 
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series.  Such  a  treatment  is  called  a  collocation  method,  but  it  is  also  referred  to  as  a 


pseudospectral  method  for  reasons  described  below.  Note  that  the  material  moduli  and 
density  already  are  sampled  at  discrete  points  in  space  by  using  an  FFT  to  perform  the 
convolutions  of  the  previous  section. 

Let  the  continuous  space  x  =  xpnp  be  discretized  into  the  positions  jpAxpnp,  with  Np 
collocation  points  evenly  spaced  by  a  distance  A xp  along  the  direction  /?.  The  wavenumbers 
become  discretized  as  -xfnP  f°r  kp  =  —Np/ 2+ 1, ...,  Np/ 2,  and  Xp  =  NpAxp.  The  spatially 
discretized  velocity  is 


N/2 

va(jpAxpnp,t)  =  Va(},t)=  J2  14(k,  t)  e'2*kplfi/Nfi  ;  jp  =  0, Np-\  (34) 

k=-N/2+l 

where  the  index  kp  =  Np/ 2  corresponds  to  both  of  the  identical  positive  and  negative  Nyquist 
frequencies  along  the  direction  f).  The  symbol  for  the  expansion  coefficient  has  been  changed 
from  V  to  V  to  distinguish  its  relationship  to  the  spatial  domain.  While  V  is  found  from 
equation  (3),  V  is  found  from  the  discrete  orthogonality  relation 


1  y;1  jihckj/N  [  1  if  j  =  nN;  n  =  0,  ±1,  ±2, ... 

N  \  0  otherwise 

so  that  the  expansion  coefficients  are  obtained  from  a  discrete  transform: 

V;(k,()=  77~  £  V.(j, t )  e-a-w/w.  ;  VN  =  N, 

Vn  j=0  T=1 


(35) 


(36) 


The  continuum  field  is  represented  by  the  N/2-degree  trigonometric  interpolant  of  the  nodal 
quantities  of  equation  (34). 

N/2 

=  J}  Va{k  (37) 

k=-N/2+l 
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Field  derivatives  in  the  discrete  space  are  defined  in  terms  of  this  continuous  field: 

jfrMx.i)  =£V„(k,0<i2**«,)'=a’W**  (38) 

k 

so  that 

=EK.(k,0(i2»W<;a*W».  (39) 

The  discrete  Fourier  expansion  coefficients  Va(k,  t)  may  be  regarded  as  approximations 
to  the  continuum  field  coefficients  V^(k,  t)  of  the  previous  section,  where  the  trapezoidal  rule 
is  used  to  evaluate  the  integral  in  the  inverse  transform. 

The  stress  field  will  be  aliased  if  the  resulting  bandwidth  of  the  convolution  sum  in  (33) 
exceeds  the  bandwidth  of  the  basis  set  used  to  synthesize  the  spatial  domain  stress  and  strain 
fields.  Let  the  spatial  strain  field  be  composed  of  a  total  of  Np  nonzero  wavenumbers  in  the 

direction  fi.  Since  we  require  that  the  stress  and  strain  fields  have  the  same  bandwidth,  the 

/  / 

indices  k  and  k  in  equation  (33)  have  the  same  range.  The  range  of  the  index  difference  kp—kp 

on  the  modulus  is  then  —Np  to  Np,  and  a  stress  field  bandlimited  to  dtNp/2  wavenumbers 
samples  the  modulus  spectrum  up  to  ±Np  wavenumbers.  If  the  bandwidth  of  the  modulus 
exceeds  ±Np,  then  the  convolution  in  (33)  will  be  aliased.  For  the  collocation  method, 
aliasing  from  the  convolution  cannot  be  avoided  if  the  bandwidths  of  the  material  structure 
and  the  wavefield  are  the  same.  The  term  pseudospectral  was  used  by  Orszag  [12]  to  describe 
such  a  method,  because  with  aliasing  the  method  is  not  a  complete  spectral  method.  Since 
our  differential  equations  are  linear  in  the  absence  of  external  forces,  however,  the  error  of 
one  wavenumber  does  not  affect  the  error  of  another  wavenumber,  and  aliasing  errors  in  the 
pseudospectral  method  are  insignificant  for  wavefields  with  most  of  their  energy  below  about 
half  the  Nyquist  sampling  frequency,  i.e.  wavelengths  less  than  about  four  grid  spaces  [15]. 
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To  obtain  a  collocation  formulation  of  the  Fourier-Sawtooth  method,  apply  the  discrete 
transform  to  equations  (23)  -(32)  to  obtain  the  following  governing  equations: 


E[  ^fn(k,0+ 

k 

(40) 

/>0)|v2Cj,()  = 

e!  ^  r„(k,t)]  ci7wk*’*iN* 

k 

(41) 

Ill 

(42) 

(43) 

The  constitutive  relation  becomes 


[  A0)  +  2/*(i)  ]  E  ^  HO. 

1 

+  A0)  [  E  +  2^V2(i„/,l)lA  ] 

/'(j)[  E(  ^  mt)  +  ^  Vk(l,0 

(44) 

(45) 

a(j)Et f^vi(i,(Kw"»/iV» 

1 

+  [AO)  +  2/*U)j  [  E  ^Vi(l.0eo""'/''»  +  3^H(j„ 

/,<)<„]  (46) 

[  AO,,/)  +  ]  E  ^ 

If 

(47) 

/.()., /)E  ‘x1  ■ 

1] 

(48) 

A(j„/)  E  TfT1  V1((„/,IKW“,/''' 

h 

(49) 

Body  forces  and/or  surface  tractions  are  applied  as  initial  conditions  on  V^(j =  0)  and 
TO0(j,  <  =  0),  and  equations  (40)  -  (49)  are  integrated  in  time  with  a  leapfrog  method 
whose  time  step  is  small  enough  to  make  numerical  dispersion  insignificant  (  e.g.,  in  1-D 
CmaxAt/Ax  ~  0.2  ).  Nyquist  errors  are  eliminated  by  using  odd-based  real-to- complex 
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FFT8  to  compute  terms  on  the  right  hand  sides  of  the  equations.  The  required  number  of 
FFT  operations  is  about  one  fourth  the  required  number  in  a  Chebychev  pseudospectral 
method  that  uses  a  fourth-order  Runge-Kutta  time  integration  scheme  and  real-to-complex 
FFTs. 

2.6.  Numerical  Tests 

The  best  test  of  the  accuracy  of  a  numerical  method’s  simulation  of  a  free  surface  con¬ 
dition  is  a  comparison  of  the  numerical  and  analytic  solutions  to  Lamb’s  problem:  An 
impulsive  source  on  the  surface  of  a  homogeneous  halfspace.  We  have  compared  analytic 
solutions  to  Lamb’s  problem  to  those  obtained  using  the  Fourier-Sawtooth  pseudospectral 
algorithm  described  in  the  previous  section.  For  the  following  comparisons,  the  P-wave  and 
S-wave  velocities  of  the  medium  are  5  km/s  and  3  km/s,  respectively,  with  a  density  of  2.5 
g/cm3,  and  the  numerical  grid  spacing  is  1  km.  For  these  values,  the  frequency  of  an  S-wave 
traveling  with  the  Nyquist  wavelength  of  2  km  is  1.5  Hz,  and  the  Rayleigh  wave  velocity  is 
about  2.74  km/s.  The  source  was  applied  as  a  delta  function  in  time  and  space  and  the  time 
series  solutions  were  lowpass  filtered  to  remove  Gibbs  truncation  effects. 

For  a  source  applied  directly  on  the  surface,  i.e.  at  x2  =  0,  considerable  energy  would 
be  transmitted  to  the  bottom  of  the  grid  because  the  basis  functions  cannot  distinguish  the 
boundary  at  x2  =  0  from  the  boundary  at  x2  =  X?.  Figure  (2)  displays  the  normalized 
kinetic  energy  density  field  of  the  Fourier-Sawtooth  solution  5.5  seconds  after  an  impulsive 
source  was  applied  at  a  depth  of  1  km.  The  field  has  been  spatially  lowpass  filtered  to 
remove  the  Gibbs  noise.  The  source  was  applied  over  four  nodes  in  the  middle  of  the  side 
of  the  grid  with  the  highest  elevation  in  the  figure.  The  P-wave,  S-wave,  and  surface  wave 
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S*w1o«h  ud  ila 
Fourier  expuaoa 


Figure  2:  Normalized  kinetic  energy  field  of  the  Fourier-Sawtooth  solution  for  an  impulsive 
source  at  a  depth  of  lifcm.  The  source  was  applied  in  the  middle  of  the  side  of  the  grid  with 
the  highest  elevation  in  the  figure. 
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phases  are  clearly  visible.  Notice  that  the  top  and  bottom  of  the  grid  are  still  coupled,  and 
surface  wave  energy  propogates  along  the  bottom  boundary  with  nearly  the  same  amplitude 
as  the  surface  waves  on  the  top  boundary.  Figures  (3)  through  (6)  compare  horizontal  and 
vertical  displacements  from  the  Fourier-Sawtooth  solution  and  a  similarity  solution  [5]  for 
the  impulsive  source  at  a  depth  of  1  km.  The  Fourier-Sawtooth  domain  has  a  uniform  grid 
spacing  of  1  km  in  each  coordinate  direction.  The  numerical  solution’s  body  wave  accuracy 
is  essentially  fhat  of  the  Fourier  method.  However,  the  numerical  solution’s  Rayleigh  wave 
is  accurate  near  the  source  but  becomes  less  accurate  as  the  Rayleigh  wave  propogates  away 
from  the  source.  The  amplitudes  of  both  the  horizontal  and  vertical  displacement  traces 
decrease  because  of  the  transmission  of  energy  through  the  boundary,  and  the  horizontal 
trace  arrives  very  early. 

The  Rayleigh  wave  solution  can  be  improved  significantly  by  refining  the  spatial  resolution 
of  the  computational  domain  in  the  vicinity  of  the  boundary.  A  spatial  domain  with  evenly 
spaced  gridpoints  is  mapped  to  a  domain  that  is  compressed  in  the  viscinity  of  the  free 
surface,  as  represented  in  Figure  (7).  Grid  mapping  in  pseudospectral  computations  was 
first  introduced  to  improve  the  accuracy  of  locating  interfaces  [7].  Tal-Ezer,  et.  al.  [13]  have 
used  a  mapping  to  increase  the  grid  spacing  near  the  boundaries  in  a  Chebychev  method 
in  order  to  allow  for  larger  time  integration  steps.  Our  mapped  spatial  domain  resembles 
the  domain  intrinsic  to  a  Chebychev  basis,  but  it  admits  a  simple  inverse  mapping  so  that 
source  and  receiver  positions  may  be  specified  easily.  Figure  (8),  when  compared  to  Figure 
(2),  indicates  that  this  mapping  has  greatly  reduced  the  amount  of  surface  wave  energy  that 
leaks  through  the  boundary.  The  leakage  is  reduced  for  body  waves  also,  although  this  is 


17 


Figure  3:  Comparison  of  horizontal  and  vertical  displacements  from  the  Fourier- Sawtooth 
solution  and  the  analytic  solution  to  Lamb’s  problem  for  a  1  km  deep  source.  The  horizontal 
and  vertical  grid  spacings  are  both  1  km,  the  source-receiver  distance  is  25  km,  and  the 
solutions  were  lowpass  filtered  with  a  corner  frequency  of  0.5  Hz. 
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Figure  4:  Same  as  Figure  (3)  except  the  source-receiver  distance  is  50  km. 
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Figure  5:  Same  as  Figure  (3)  except  the  source- receiver  distance  is  100  km 
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Figure  6:  Same  as  Figure  (3)  except  the  source-receiver  distance  is  200  km. 
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not  apparent  given  the  amplitude  scale  of  the  figure.  Notice  that  some  energy  in  the  grid 
is  still  trapped  in  the  source  region  at  the  time  of  the  snapshot.  Because  the  source  was 
applied  locally  in  a  region  with  relatively  small  grid  spacings,  it  generated  high-  frequency 
energy  that  cannot  be  supported  by  the  larger  grid  spacings  at  depth.  These  frequencies 
are  reflected  and  trapped  within  the  finer  grid.  The  simulation  used  to  generate  Figure 
(8)  was  identical  to  the  one  for  Figure  (2)  except  for  the  domain  mapping.  The  initially 
uniform  domain  with  a  grid  spacing  of  1.0  km  was  mapped  to  a  domain  with  a  minimum 
grid  spacing  of  0.20  km  at  the  boundaries  and  a  maximum  spacing  of  1.89  km  at  the  center 
of  the  domain.  Figures  (9)  through  (12)  compare  displacements  from  the  Fourier-Sawtooth 
and  analytic  solutions  for  a  1  km  deep  source  and  this  domain  mapping  in  the  vertical 
direction  of  the  numerical  method.  The  horizontal  grid  spacing  was  kept  uniform  at  1  km. 
The  horizontal  displacement  solution  has  improved  significantly,  and  both  horizontal  and 
vertical  displacement  solutions  decay  very  ;  'ttle  with  distance  from  the  source.  Both  traces’ 
arrival  times  are  in  error  by  slightly  more  than  one  percent.  This  is  to  be  expected  for  an 
approximation  method  based  on  a  variational  principle.  The  approximate  eigenvalues  are 
slightly  higher  than  the  exact  ones.  In  Figure  (12),  the  amplitude  mismatch  after  74  seconds 
is  due  almost  entirely  to  interference  from  a  body  wave  incident  from  depth,  as  no  absorbing 
boundaries  were  applied. 

We  have  compared  the  Fourier-Sawtooth  solution  to  a  normal  mode  solution  [9]  for  a  1 
km  deep  explosive  source  in  the  structure  shown  in  Figure  (13).  A  high-velocity  cap  layer 
was  included  in  the  structure  at  a  depth  of  150  km  for  the  normal  mode  simulation.  Ab¬ 
sorbing  boundary  conditions  were  applied  to  the  Fourier-Sawtooth  algorithm  by  the  method 
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suggested  by  Cerjan,  et.  ai.  [2].  The  velocity  and  stress  fields  are  attenuated  within  a  zone 
of  grid  points  near  the  boundaries,  and  the  amount  of  attenuation  increases  as  the  wavefield 
approaches  the  boundaries.  The  results  were  lowpass  filtered  with  a  corner  frequency  of 
0.5  Hz.  The  time  series  from  the  two  methods  at  a  source-receiver  distance  of  205  km  are 
overlayed  in  Figure  (14)  to ‘compare  the  surface  wave  dispersion,  and  record  sections  from 
the  two  algorithms  are  shown  in  Figure  (15).  Because  the  normal  mode  algorithm  has  cylin¬ 
drical  symmetry  and  the  Fourier- Sawtooth  method  is  cartesian,  the  relative  body  wave  to 
surface  wave  amplitudes  do  not  match,  but  the  surface  waves  display  very  similar  dispersion 
characteristics  when  scaled  for  comparison. 
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Figure  9:  Comparison  of  horizontal  and  vertical  displacements  from  the  Fourier-Sawtooth 
solution  and  the  analytic  solution  to  Lamb’s  problem  for  a  1  km  deep  source.  The  horizontal 
grid  spacing  is  constant  at  1  km  but  the  vertical  spacing  is  mapped  to  refine  the  spatial 
resolution  in  the  vicinity  of  the  boundary.  The  source-receiver  distance  is  25  km,  and  the 
solutions  were  lowpass  filtered  with  a  corner  frequency  of  0.5  Hz. 
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Figure  12:  Same  as  Figure  (9)  except  the  source- receiver  distance  is  200  km.  The  amplitude 
mismatch  after  74  seconds  is  due  almost  entirely  to  interference  from  a  body  wave  incident 
from  depth. 
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Figure  (13).  Record  sections  of  the  Fourier-Sawtooth  and  Normal  Mode  solutions  for  the 
structue  shown  in  Figure  (13). 


2.7.  Numerical  Modeling  of  Seismic  Sources  and  Wave  Propagation  in  Complex  Media 

We  have  recently  developed  several  numerical  modeling  programs  that  are  well  suited  to 
both  non-linear  source  modeling  problems  and  wave  propagation  in  elastic-anelastic  or  plastic 
media.  These  finite  difference  and  finite  element  type  programs  are  all  operational  on  our 
Stardent  computer  system  which  allows  us  to  perform  lengthy  computations  in  2  and  3  dimen¬ 
sional  grid  systems.  The  graphics  features  of  the  Stardent  are  particularly  valuable  in  display¬ 
ing  and  understanding  the  results. 

The  particular  numerical  programs  we  have  available  for  this  study  are  of  two  classes; 
those  that  can  be  called  the  "standard”  finite  difference  or  finite  element  programs,  and  pro¬ 
grams  based  on  spectral  or  pseudospectral  methods  {eg.  Carcione  et.  al.  1992,  Kosloff  et.  al. 
1990,  Witte  and  Richards,  1990,  Canuto  et.  al.  1988,  Orszag  1971).  The  "spectral"  programs 
have  been  developed  recently  {eg.  Orrey  and  Archambeau,  1993)  using  generalizations  of  the 
earlier  pseudospectral  methods,  as  well  as  refinements  of  numerical  procedures,  to  produce  a 
fast,  memory  economic  and  accurate  computational  method  for  the  simulation  of  seismic  wave 
propagation  in  2  and  3-D  anelastic  media. 

Tests  of  the  "full  spectral"  numerical  method  (wherein  the  spatial  dependence  of  the  field 
variables,  as  well  as  the  medium  structure  properties,  are  represented  in  terms  of  Fourier  basis 
sets)  and  the  psuedospectral  method  (wherein  only  spatial  derivatives  of  field  variables  are 
computed  by  FFT  methods)  have  been  made  to  check  the  accuracy  of  the  methods.  Our 
approach  has  been  to  compare  these  numerical  results  to  those  obtained  by  modal  synthesis 
methods  {eg.  Harvey,  1981)  for  laterally  uniform  models.  In  this  regard  Figure  (16)  shows  a 
comparison  of  synthetics  produced  by  the  analytically  based  modal  synthesis  method  and 
those  by  the  pseudospectral  method. 
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The  structure  used  is  shown  in  Figure  (17),  and  is  applicable  to  the  region  of  the  nuclear  test 
site  in  Eastern  Kazakhstan.  Similar  results  are  obtained  using  the  "full  spectral"  method 
described.  These  tests  show  very  good  accuracy  for  these  methods  based  on  computational 
comparisons  to  modal  (and  reflectivity)  methods.  Consequently,  we  feel  confident  that  these 
(new)  numerical  methods  will  also  be  accurate  for  wave  propagation  studies  in  the  laterally 
variable  2  and  3-D  models  .of  interest. 

In  this  regard  Figure  (18)  illustrates  results  from  full  2-D  tests  of  these  modeling  pro¬ 
grams.  In  this  example  we  have  examined  some  effects  of  near  source  surface  topography 
and  fine  scale  layering  on  the  seismic  wave  field  generated  by  a  buried  explosion  (depth  300 
meters).  The  sequence  of  insets  in  Figures  (18a)  and  (18b)  shows  die  spatial  evolution  of  the 
wave  field  with  time  following  die  detonation  of  the  explosion,  with  the  effects  of  topographic 


Figure  17.  Layered  structure  used  in  the  modeling  comparisons  of  Figure  1 . 
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(a.)  Vertical  Velocity  Field  32  sec.  after  detonation 


(b.)  Vertical  Velocity  Field;  38  sec.  after  detonation 


(c.)  Vertical  Velocity  Field;  .55  sec  after  detonation 


Figure  18-a.  Contours  of  the  vertical  particle  velocity  from  an  explosion  300  meters  below  the  free  surface 
with  histh  relief  surface  topography  in  the  vicinity  of  the  source  area  The  dimensions  of  the 
err-,.; -section  shown  is  9  km.  wide  by  3. 45  km.  deep.  Maximum  surface  elevation  is  300  meters 
Red  denotes  zones  of  a  downward  vertical  velocity  component.  Mue  is  ar,  upward  velocity 
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relief  on  the  surface  reflections  following  the  direct  P  wave,  for  example,  producing  rather 
large  variations  in  the  P  wave  amplitudes  as  a  function  of  emergence  angle  from  the  source 
area.  These  effects  can  also  be  seen  in  the  time  series  recorded  at  the  free  surface  at  greater 
distance  from  the  source  zone.  In  this  regard.  Figure  (19)  shows  recorded  particle  velocity  at 
6  km.  from  the  source,  under  conditions  in  which  the  rough,  near-source  topography  is  either 
present  (b.  and  d.)  or  absent  (a.  and  c.).  The  effects  of  shallow  fine  scale  layering  are  also 
indicated  by  these  examples. 


(t.)  Vertical  Vckxajr  ii»e  *na  acar  llie  ««f«cc  cf  «n  cUflic 
hommgemtoms  half  gMC  *******  aarface 


(t )  Vertical  Velocity  wme  xcncs  near  the  surface  of  a  (5)  lay<*td 
half  space  without  surface  topography 


(bj  Vertical  Velocity  lane  aeries  near  the  surface  of  an  elastic 
komofauoms  Uf  space  with  rough  strfsoc  topogra¬ 

phy  near  the  aouroc  area. 


(d.)  Vertical  Velocity  bme  series  near  the  surface  of  a  (5)  loytrfd 
half  space  with  rough  surface  topography  near  the  source  area 


Figure  19.  The  effects  of  near  source  topography  and  fine  scale  layering  on  the  wave  field  from  a 
buried  explosion  source  (300  m  depth).  Observations  are  at  6  km  from  the  source  on  the  free  surface. 
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Clearly  the  examples  show  that  complexities  in  the  P  and  Pn  coda,  of  the  type  observed, 
can  be  produced  by  topographic  features  and/or  near  source  fine  scale  layering.  In  future  stu¬ 
dies  we  will  systematically  examine  these  effects,  using  realistic  topography  and  structural 
layering  and  compare  these  results  to  observational  data.  In  addition,  using  extended  versions 
of  the  numerical  modeling  programs,  whose  out-put  is  illustrated  here,  we  will  include  failure 
phenomena  in  the  near  source  region  to  produce  spallation,  in  order  to  evaluate  its  contribu¬ 
tion  to  the  radiated  wave  field. 

The  effects  of  strong  and  moderate  lateral  variations  in  near  source  shallow  structure  can, 
of  course,  also  be  included  and  can  be  systematically  studied  along  with  all  the  other  medium 
characteristics  producing  strong  perturbations  in  the  wave  field.  Through  this  approach  we 
hope  to  be  able  to  "sort  out"  the  quantitative  nature  of  each  of  the  effects  on  the  directly  radi¬ 
ated  seismic  field  from  explosions  and  be  able  to  evaluate  their  total  impact  on  source  depth 
estimates,  yield  estimation  and  discrimination.  As  noted,  the  effects  of  strong  "random" 
fluctuations  in  material  properties  at  upper  and  mid-crustal  depths  is  also  of  importance  in 
producing  scattering  that  can  enhance  Pg  and  Lg  excitation,  while  reducing  Rg  by  scattering 
losses  with  distance.  Figure  (20)  shows  our  preliminary  modeling  tests  using  numerical  simu¬ 
lations  in  2-D  models  with  randomized  intrinsic  velocities.  Clearly,  the  complex  forms  of 
observed  seismic  data  have  a  strong  resemblance  to  the  synthetics  obtained  in  a  randomized 
earth  model.  Of  course  our  approach  will  be  to  systematically  investigate  this  kind  of  scatter¬ 
ing  effect,  particularly  as  a  loss  mechanism  for  Rg  and  as  an  excitation  mechanism  for  Lg  and 
Pg,  while  also  carefully  evaluating  the  "trade-off  of  anelastic  attenuation  versus  this  scattering 
process  for  Rg,  as  well  as  the  other  possible/probable  mechanisms  of  Pg  and  Lg  excitation  due 
to  major  lateral  variations  in  structure. 
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Figure  20.  Effects  of  vertical  elastic  parameter  randomization,  insets  a.  and  b.  are  vertical  and  radial 
velocity  components  without  radomizatioo.  Insets  c  and  d.  are  for  the  same  structure  with  superim¬ 
posed  randomization  of  25%  at  the  surface  to  10%  at  the  base  of  the  crust  for  the  elastic  velocity  vari¬ 
ation  with  depth. 
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In  regard  to  the  latter  effects  involving  major  structural  discontinuities.  Figure  (21)  shows 
examples  of  the  effects  of  a  major  shallow  structural  transition,  in  this  case  a  deep,  low  velo¬ 
city,  sedimentary  basin  within  a  higher  velocity,  largely  igneous,  crust.  As  illustrated,  features 
of  this  type  strongly  affect  high  frequency  surface  waves,  like  Rp  and  produce  mode  coupling 
effects  that  give  rise  to  enhanced,  and  complicated,  Pg  and  Lg  excitation.  We  therefore  will 
study  such  strong  structure  induced  effects  on  a  systematic  basis  to  try  to  quantify  the  varia¬ 
bility  of  discriminatory  signals,  like  Pg  and  Lg,  in  order  to  assess  what  characteristics  of  these 
signals  nevertheless  remain  robust  and  could  be  used  for  reliable  discrimination  under  a 
variety  of  structural  conditions. 

Thus,  with  our  present  analytical  modeling  capabilities  for  both  sources  and  wave  propa¬ 
gation,  in  both  vertically  and  laterally  varying  media,  we  expect  to  provide  a  means  of 
predicting  and  understanding  very  complex  source  and  near  -  source  behavior  (eg.  multiple 
explosions  in  fractured  media)  as  well  as  complexities  in  wave  propagation  effects. 
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3.  Modeling  Atmospheric  Wave  Fields  and  Ionospheric  Electron  Density  Variations  Due 
to  Near  Surface  Seismic  Sources 

3.1.  Modeling  of  Atmospheric  and  Ionospheric  Gravity  Waves 

Because  of  the  exponential  decrease  of  atmospheric  density  with  height,  buoyant  pulsed 
gravity  waves  generated  by  surface  or  subsurface  seismic  sources  can  be  of  appreciable  ampli¬ 
tude  through  out  the  atmosphere.  Furthermore,  above  100km  in  height,  these  flow  transients 
affect  the  ionospheric  EM  fields  through  changes  in  the  distribution  of  the  charged  particles. 
Accurate  modeling  estimates  of  seismically  generated  gravity  waves  and  their  effects  in  the 
ionosphere  have  not,  to  our  knowledge,  been  investigated. 

The  basic  equations  governing  motions  of  the  neutral  atmosphere  are  the  conservation 
laws  of  mass,  momentum  and  energy  together  with  the  ideal  gas  equation  of  state.  The 
specific  nonlinear  continuum  equations  incorporate  advective  terms  as  well  as  the  gravitational 
field,  gas  compressibility,  viscosity  effects  and  thermal  conductivity.  For  electron  motions  in 
the  ionosphere  a  first-order  continuity  equation  is  used  which  assumes  that  electrons  move 
with  the  neutral  atmosphere. 

In  our  work  to  date,  the  set  of  partial  differential  equations  for  the  atmosphere  are  con¬ 
verted  to  a  corresponding  set  of  finite  difference  equations  in  order  to  effect  numerical 
integration  in  time  and  space.  The  non-linear  terms  are  treated  non-locally  on  the  lattice  for 
stability,  effectively  controlling,  internally,  the  instabilities.  In  addition,  random  velocities  and 
pressures  are  attributed  to  the  inherent  fine  scale  turbulence  in  the  atmosphere  and  are  incor¬ 
porated  in  the  modeling  as  are  mean  drift  particle  velocities.  In  particular,  in  order  to  account 
for  the  inherent  turbulence  in  the  atmosphere,  the  flow  variables  at  a  point  are  decomposed 
into  a  mean  flow,  governing  winds,  and  a  perturbed  flow  that  incorporates  turbulence.  A  new 
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approach,  designed  to  include  turbulence,  has  been  developed  using  random  perturbations 
obtained  from  a  random  number  generator  which  are  input  directly  into  the  finite  difference 
equations.  Turbulence  is  also  produced  by  a  random  distribution  of  temperature  at  the  surface 
which  produces  thermal  structures  with  upward  and  downward  flows.  Horizontal  winds, 
impacting  on  a  variable  and  random  topography,  also  produce  upward  and  downward  motions 
which  have  a  random  stochastic  character. 

The  set  of  finite  difference  equations  are  numerically  integrated  in  time  and  space. 
Upwind  differencing  is  used  for  first  order  spatial  gradients  with  the  advection  velocity  terms 
acting  at  the  upwind  point.  However,  if  the  velocity  operates  on  its  own  velocity  gradient, 
such  non-linear  terms  are  treated  non-locally  on  the  lattice  for  stability,  effectively  controlling 
internally  any  unstable  growth. 

There  are  at  least  three  types  of  boundary  important  in  this  modeling.  The  air-ground 
surface  is  topographically  complex  with  a  turbulent  boundary  layer  of  the  order  of  a  few 
meters  at  the  interface.  At  this  boundary,  vertical  velocities  are  random  both  in  time  and  at 
spatial  locations.  Because  of  the  presence  of  the  lower  boundary  layer  above  a  complex 
topography,  horizontal  velocities  are  not  taken  as  zero  but  incorporate  winds  and  turbulence 
effects.  The  top  atmospheric  boundary  is  open  with  decreasing  density.  The  topmost  boun¬ 
dary  should  mimic  the  conditions  for  an  open  atmosphere  with  specific  considerations  for 
buoyancy  and  field  gradients.  We  have  examined  various  options  including  fixing  velocities 
and  densities  and  their  gradients.  However,  we  have  adopted  the  general  open  flow  boundary 
such  as  we  also  used  for  the  artificial  side  boundaries.  The  side  boundaries  are  artificial,  due 
to  grid  restrictions,  and  must  mimic  open  boundaries  that  allow  free  flow  in  either  direction. 
We  have  adopted  the  more  usual  approach  wherein  the  dependent  variables  are  constrained  to 
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stay  constant  at  these  open  boundaries. 


3.2.  Conservation  Laws 


The  continuity  equations  are  based  on  the  values  of  the  fields  at  particular  points  in 
space  and  time.  Conservation  of  mass  is  expressed  as, 


(1.) 


where  p  is  density  and  Uj  is  velocity  in  the  Xj  direction.  Conservation  of  momentum  is  simi¬ 
larly  expressed  as: 


dU)  chi; 


dt  +Uj  dx 


‘JJ 


(2.) 


where  Xj  are  external  forces,  Py  is  the  generalized  stress  such  that: 


with  ey  =  1/2 


duj  du: 

-  +  — L 

dxj  dx; 


Py  -  -V  ’  &ij  +  •  esj  -  2/3p  •  6jj  •  (3.) 

the  strain  rate  and  |i  the  viscosity.  Conservation  of  energy  is: 


dT 

dxj 


du: 

-P '  aT  +  <i> 

OX: 


(4.) 


where  T  is  temperature,  Cy  is  specific  heat  at  constant  volume,  K  is  thermal  conductivity,  and 
<t>  is  the  viscous  dissipation.  Here: 

*-2lie}-2 /3-lKe*)2 

The  equation  of  state  for  the  atmospheric  gas  is  taken  to  be  ideal,  i.e. 


P  «  ^5.  .  p  .  t 

m 


(5.) 


where  kB  is  Boltzmann’s  constant  and  m  is  the  mean  molecular  weight. 
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3.3.  Normalized  Equations 


In  the  fundamental  equations,  (1)  through  (5),  the  dependent  and  independent  variables 
can  be  normalized  with  respect  to  typical  values.  For  an  ambient  atmosphere,  with  exponen¬ 
tial  decay  of  density  with  height,  distances  are  normalized  through  the  scale  height,  H,  which, 
at  the  surface,  is  approximately  8400  metres.  Velocities  are  normalized  with  respect  to  c,,  the 
sound  velocity  of  air  at  the  earth’s  surface.  Similarly  density,  pressure  and  temperature  are 
normalized  to  surface  values,  and  the  independent  variable,  time  t,  is  normalized  by  (H/c„). 
Thus  for  the  continuity  of  mass,  we  get,  as  before, 

f  +  A(pUj).o  (6.) 

where  the  new  variables  are  now  normalized  and  given  the  same  symbol  as  the  original  vari¬ 
ables.  Incorporating  gravity  as  the  external  force,  the  momentum  conservation  equation,  (2), 
becomes 


p~dt  +  pu> '  — G* -g(z)p  ^"A2  ^  +  (7) 

where  G„  *  (gJH/c,2)  is  a  measure  of  the  ratio  of  potential  energy  to  thermal  energy. 

A2  =  PsKPgC*2)  is  a  measure  of  the  ratio  of  stress  energy  to  thermal  energy. 

A4  *  Ps/lpj  Cg)  is  the  ratio  of  viscous  to  thermal  energy. 

'P  is  the  normalized  viscosity  drag  force, 


8uj  duj 
dxj  dxj 


duk 


where  p  is  normalized  to  p,,  the  viscosity  at  the  surface.  In  the  atmosphere,  p  is  usually 


taken  to  be  constant  for  the  molecular  viscosity.  In  the  case  of  conservation  of  energy,  the 
normalized  equation  is 
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dT  dT  .  d*T  A  dUi  .  . 
3t  Uj'  dXj  “A®'  dx,2  ~A*'  a*  +  As* 


where 


As-H. 


dXj 


^  Ca  ^  .  \ 
m-Cv  (p,H)  A2 


(8.) 


m-Cy 


A*K 


(P.C.H) 

where  K  is,  as  usual,  taken  constant  for  the  atmosphere.  The  equation  of  state,  on  normaliza¬ 
tion,  is 


p  «  p  •  T  •  /m(z)  (9.) 

kn 

where  the  ideal  equation  of  state  at  the  surface  is  p, ■  — p,T,  with  m,  the  surface  value  of 

no, 

mean  molecular  weight  (29.0)  and  m(z)  the  height-dependent  normalized  value. 

3.4.  Ionospheric  Motions 

The  basic  conservation  law  of  charged  particles,  assuming  no  creation  or  annihilation,  is: 

(10.) 

(ft  OXj 

where  Na  is  the  number  of  particles  of  type  a  and  Uja  is  their  velocity  in  the  j’th  direction. 
The  initial  concentration  of  the  charged  particles  is  taken  to  be  time-independent,  with  only  a 
vertical  functional  dependence,  N(z),  where  the  subscript  a  has  been  dropped  for  the  type  of 
particle.  Assuming  only  small  changes  in  this  concentration,  the  dependence  can  be  found 
from  integrating  eqn.  (10)  over  the  range  to  to  the  present.  To  zeroth  order,  this  concentration 
change  becomes: 
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6N(z,t)  -  -  •  |  utfW  -  N(z)  •  |  ^pdt'  (11.) 

The  firs t  term  in  (11)  is  the  concentration  change  due  to  the  displacement  of  the  ionospheric 
layer,  while  the  second  term  arises  as  a  result  of  compression  or  rarefaction  and  is  the 
predominant  term  when  dealing  with  processes  involving  characteristic  dimensions  smaller 
than  the  width  of  the  layer.  The  velocity  of  the  charged  particle  is  usually  assumed  to  be 
identical  with  that  of  the  neutral  gas  to  zeroth  order  and  this  is  the  velocity  that  is  used  in  the 
finite  difference  calculations  for  electron  density  changes. 

The  initial  concentration  of  electrons  is  taken  to  be  that  of  a  Chapman  distribution  which 
has  a  maximum  density  at  345  km  height  and  decreases  rapidly  below  about  90  km  with  the 
functional  dependence  on  height  defined  by: 

N(z)  =  Nc  •  exp  (-|(K-e^))  (12.) 

Where  *  (z-hc)/H.  h«.  =  345  km,  H  =  65  km  and  Nc  is  the  normalizing  value. 

3.5.  Finite  Difference  Scheme 

The  set  of  non-linear  partial  differential  equations  are  converted  to  a  corresponding  set  of 
finite  difference  equations  for  explicit  computer  integration  in  time  and  space.  Upwind 
differencing  is  used  for  first  order  spatial  gradients,  with  the  advection  velocity  terms  acting  at 
the  upwind  point.  However,  if  the  velocity  operates  on  its  own  velocity  gradient,  such  non¬ 
linear  terms  are  treated  non-locally  on  the  lattice  for  stability,  effectively  controlling  internally 
any  unstable  velocity  growth. 

The  updated  variable  is  projected  not  from  just  the  old  dependent  variable,  a  process  that 
is  inherently  unstable,  but  from  a  distributed  smoothed  average  of  the  variable  at  locations 
surrounding  the  specific  spatial  location.  Such  a  smoothing  method  brings  stability  to  the 
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differencing  scheme.  However,  the  attendant  numerical  diffusion  is  minimized  by  not 
smoothing  the  density  variable,  which  has  only  a  small  effect  on  stability.  This  approach  also 
helps  in  stabilizing  the  integration  at  grid  corners  and  boundaries.  The  second  order  deriva¬ 
tives  in  the  viscosity  and  thermal  conductivity  terms  are  modeled  by  finite  differences  taken  at 
the  surrounding  spatial  locations. 

In  the  explicit  integration  scheme,  the  updated  flow  velocities,  temperature  and  density 
are  obtained  via  their  continuity  equations  while  pressure  is  obtained  from  insertion  of  the 
updated  density  and  temperature  into  the  ideal  gas  equation. 

3.6.  Boundary  Conditions 

There  are  at  least  three  types  of  boundaries  important  to  the  modeling  of  fluid  flows. 
The  air-ground  surface  is  topographically  complex  with  a  turbulent  boundary  layer  attached. 
The  top  atmospheric  boundary  is  open  to  space  with  decreasing  density.  The  side  boundaries 
are  artificial,  due  to  grid  restrictions,  and  must  mimic  open  boundaries  that  allow  free  flow  in 
either  direction. 

At  the  bottom  boundary  vertical  velocity  functions  are  input  as  sources  of  m  omenta  at 
spatial  locations.  Otherwise,  as  usual,  vertical  velocities  are  taken  to  be  zero  at  the  bottom. 
Because  of  the  presence  of  the  lower  boundary  layer  above  a  complex  topography,  horizontal 
velocities  are  not  taken  as  zero  but,  instead,  constant  velocity  and  density  gradients  are 
assumed  in  the  vertical  direction.  For  subsurface  sources  only  momentum  inputs  are  con¬ 
sidered  at  this  bottom  boundary.  For  sources  at  or  above  the  surface,  both  momentum  and 
pressure  conditions  must  be  applied,  just  as  in  atmospheric  sources. 

The  top-most  boundary  should  mimic  the  conditions  for  an  open  atmosphere  with 
specific  considerations  for  buoyancy  and  field  gradients.  We  have  examined  various  options 
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including  fixing  velocities  and  densities  and  their  gradients.  However,  we  have  adopted  the 
general  open  flow  boundary,  much  as  we  use  for  the  artificial  side  boundaries.  As  usual,  in 
order  to  preserve  conservation  relations,  all  normal  gradients  are  set  to  zero  at  these  open 
boundaries.  However,  this  would  not  permit  heat  flow  through  the  boundary.  Therefore  the 
second-order  normal  derivative  of  temperature  is  made  constant. 

3.7.  Turbulence  Effects 

In  order  to  account  for  the  inherent  turbulence  in  the  atmosphere  below  the  thermopause, 
the  flow  variables  at  a  point  can  be  decomposed  into  a  mean  flow  and  a  perturbed  turbulent 
flow.  In  the  momentum  equation,  additional  components  are  thus  obtained  for  the  generalized 
stress.  These  are  mainly  interaction  terms  between  the  mean  and  perturbed  densities  and 
velocities,  termed  the  Reynold’s  stresses,  which  represent  the  interaction  of  the  mean  flow 
with  the  background  turbulence.  These  extra  stresses  have  been  approximated  by  various 
phenomenological  approaches.  Boussinesq  introduced  the  concept  of  eddy  viscosities  in  order 
to  use  the  Newtonian  equations  with  the  usual  but  much  larger  viscosity  term.  Our  models 
evaluate  the  efficacy  of  this  method  using  an  eddy  thermal  conductivity.  We  also  attempt  to 
evaluate  different  forms  of  these  Reynold’s  stresses  through  algorithmic  modeling  of  various 
drag  forces  that  mimic  the  effect  of  these  interaction  terms.  It  is  found  that  even  small  drag 
forces  of  a  particular  type  can  alter  the  flows  and  their  temporal  dependence.  An  alternative 
approach  to  turbulence  is  developed  in  the  use  of  random  perturbations,  obtained  from  a  ran¬ 
dom  number  generator,  and  input  directly  into  the  finite  difference  equations. 

3.8.  Modeling  Results 

Explosive  sources  at  and  below  the  ground  are  simulated  by  computations  in  a  solid  con¬ 
tinuum  and  their  resultant  effects  on  the  atmosphere,  at  the  solid-air  interface,  are  integrated 
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upward  and  outward.  Various  velocity  sources  are  used  at  the  lower  boundary  with  differing 
time,  amplitude  and  radial  dependences.  The  standard  input  is  a  source,  comprised  of  the  first 
differential  of  a  gaussian  in  time,  that  approximates  the  initial  pulse  from  an  underground 
explosion.  Cartesian  coordinates  are  used  to  model  the  3-dimensional  system,  with  the  source 
at  the  center  of  the  bottom  plane. 

Figure  (22)  shows  results  of  modeling  the  low  frequency  gravity  waves  (or  acoustic- 
gravity  coupled  waves)  in  the  atmosphere  and  ionosphere  due  to  a  contained  underground 
explosion  just  beneath  the  ground  surface.  The  predicted  electron  density  fluctuations,  using 
equation  (11.),  are  also  shown  in  the  upper  left  panel.  These  results  show  that  the  maximum 
of  the  disturbance  is  well  behind  the  acoustic  wave  at  the  front  which,  for  the  instant  of  time 
shown  here,  is  well  above  the  400  km.  altitude  and  is  also  small  relative  to  the  amplitude  of 
the  fields  indicated  for  this  gravity  controlled  disturbance.  Figures  (23)  and  (24)  show  the 
disturbance  at  early  and  late  times,  whereas  Figure  (22)  shows  an  intermediate  time  relative  to 
these  latter  two  "snap-shots"  of  the  fields.  Note  that  the  late  time  disturbance  is  quite  com¬ 
plex,  with  the  velocity  and  temperature  fields  showing  an  oscillatory  character  as  functions  of 
altitude.  (This  is  also  true  of  the  electron  density  but  is  not  obvious  from  the  color  coded 
amplitude  scale  used  for  the  figures.) 

In  Figure  (23),  which  shows  the  disturbance  at  an  early  time,  the  true  acoustic  wave 
front  is  somewhat  above  300  km,  well  above  the  large  amplitude  disturbance  shown  here,  and 
as  an  amplitude  lower  than  the  first  level  of  the  coarse  color  table  used  and  does  not  show  up 
above  the  background  color  in  the  figure.  In  any  case,  the  pure  acoustic  wave  is  small  rela¬ 
tive  to  the  lower  frequency,  gravity  controlled,  disturbance  shown  here.  Consequently,  it  is 
this  strongly  nonlinear  disturbance  that  appears  to  be  the  best  target  for  detection  by 
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ATMOSPHERIC  AND  IONOSPHERIC  FIELDS  DUE  TO  GROUND  MOTIONS 


ELECTRON  DENSITY  PRESSURE  TEMPERATURE 


VERTICAL  VELOCITY  HORIZONTAL  VELOCITY  INWARD  VELOCITY 


Figure  22.  CROSS-SECTIONS,  height  400km  width  800km,  of  spatial  distributions  of 
changes  in  atmospheric  and  ionospheric  fields  due  to  ground  motions  produced  by  a  subsur¬ 
face  explosion.  Acoustic-gravity  waves  propagate  upward  and  outward  in  the  ARDC  standard 
atmosphere.  Electrom  density  changes  in  a  chapman  model  are  due  to  the  electrons  following 
the  motions  of  the  predominant  neutals.  This  pattern  depicted  occurs  at  a  time  of  about  10 
minutes  after  ground  movement  with  onset  of  wave  flipping  at  the  100  km  region.  Blue 
colors  denote  negative  values  with  red  denoting  positive  values. 
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EARLY  TIMES 


ELECTRON  DENSITY  PRESSURE  TEMPERATURE 


VERTICAL  VELOCITY  HORIZONTAL  VELOCITY  INWARD  VELOCITY 


Figure  23.  CROSS-SECTIONS,  height  400km,  width  800km,  of  spatial  distributions  of 
changes  in  atmospheric  and  ionospheric  fields  due  to  ground  motions  produced  by  a  synthetic 
subsurface  explosion. 


LATE  TIMES 


ELECTRON  DENSITY  PRESSURE  TEMPERATURE 


VERTICAL  VELOCITY  HORIZONTAL  VELOCITY  INWARD  VELOCITY 


Figure  24  CROSS-SECTIONS,  height  400km,  width  800km,  of  spatial  distributions  of 
changes  in  atmospheric  and  ionospheric  fields  due  to  ground  motions  produced  by  a  synthetic 
subsurface  explosion. 
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electromagnetic  sounding  methods  since  the  associated  electron  density  fluctuations  are  many 
time  larger  than  those  from  the  acoustic  (or  "shock")  front. 

A  more  detailed  representation  of  the  field  variables  (velocity  and  pressure)  are  shown  in 
Figure  (25).  Here  the  fluctuating  nature  of  the  predicted  disturbance  is  evident.  The  general 
predicted  form  of  this  disturbance,  which  is  reflected  in  electron  density  variations,  should 
allow  this  nonlinear  wave  to  be  easily  detected  by  E-M  sounding  methods. 

The  results  of  the  atmospheric  modeling  of  the  gravity  wave  effects  from  a  surface 
explosion  can  be  summarized  as  follows: 

(1.)  A  time  dependent  transient  pulse  propagates  upward  with  increasing  amplitude  rela¬ 
tive  to  the  ambient  pressure.  This  produces  asymmetric  flows  which  control  the  flow 
development  and  the  upward  propagation  of  the  transient.  The  initial  positive  density 
pulse  is  propagated  upward  more  slowly  than  the  following  negative  density  pulse  which 
has  increased  buoyancy.  This  initiates  a  sequence  of  circulation  patterns  that  develops 
through  what  appears  to  be  asymmetric  triangular  modes  across  the  horizontal  cross- 
section.  The  circulation  patterns  for  the  phenomena  are  characterized  by  upward  central 
motions  of  the  lighter  matter,  which,  at  the  neutral  buoyancy  level,  push  outward  to  the 
side.  The  centroid  of  the  transient  pulse  initially  moves  upward  rapidly,  but  slows  down 
to  the  group  velocity  speed  of  sound  in  the  atmosphere.  The  advected  air  mass  tries  to 
remain  in  its  horizontal  stratification  in  order  to  minimize  changes  in  its  gravitational 
potential.  However,  it  appears  that  energy  and  momenta  are  transported  through  travel¬ 
ing  waves  in  the  circulation  pattern.  Similar  effects  have  been  observed  in  the  real  atmo¬ 
sphere  when  thermals  propagate  upward  from  the  Earth’s  surface  with  similar  circulation 
patterns. 
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(2.)  After  a  model  dependent  characteristic  time,  a  bifurcation  of  the  flow  occurs  with 
the  eventual  reversal  of  the  velocity  directions.  The  bifurcation  phenomena  occurs  in 
this  model,  every  100  seconds,  so  that  it  has  a  period  of  just  over  3  minutes.  A  drag 
force  is  input  in  order  to  model  the  effect  of  the  inherent  background  turbulence  of  the 
atmosphere.  A  drag  force,  which  removes  2%  of  the  component  velocities  at  each  com¬ 
putational  grid  point  at  each  time  step,  removes  the  periodic  bifurcation  and  a  standing 
wave  is  formed  in  the  atmosphere  with  constant  field  patterns.  However,  with  a  1% 
removal  rate,  die  patterns  are  periodic  with  similar  bifurcations  as  in  the  zero  drag  case. 
Because  existing  atmospheric  turbulence  acts  on  the  transient  gravity  wave  as  a  perturba¬ 
tion,  we  have  also  modeled  its  effect  by  imposing  a  random  component  on  each  field  at 
each  time  step  and  grid  point.  The  usual  bifurcations  are  obtained  but  with  differing  pat¬ 
terns  from  the  zero  turbulence  case.  However,  the  appearance  of  the  pressure  and  den¬ 
sity  fields  is  more  realistic  due  to  added  diffusion  and  random  components. 

As  the  transient  pulse  moves  upward  in  the  atmosphere,  it  magnifies  in  amplitude  relative 
to  the  exponentially  decreasing  ambient  pressure.  Thus,  the  level  at  which  a  specific  pressure 
is  located  will  oscillate  as  the  transient  pressure  pulse  moves  through.  To  the  first  order,  the 
electrons  in  the  ionosphere  are  assumed  to  move  with  the  flow  of  the  dominant  neutrals. 
Thus  the  change  in  the  electron  density  can  be  calculated  from  a  conservation  law,  whose 
integration  in  time  gives  the  total  electron  density  variation.  The  ambient  electron  density  is 
approximated  by  the  Chapman  function  which  has  a  maximum  electron  density  at  350  km  and 
effectively  zero  electron  density  below  about  90km.  For  reasonable  velocity  sources  at  the 
ground  surface,  we  find  that  changes  in  electron  density  from  100  km  and  upward  are  at  least 
of  the  same  order  as  those  observed  by  EM  experiments  conducted  over  surface  and  subsur¬ 
face  explosions,  and  generally  for  the  gravity  wave  disturbance  it  appears  to  be  larger.  In  this 
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regard.  Figure  (26)  shows  an  example  of  the  predicted  gravity  wave  induced  fluctuations  in 
temperature  and  electron  density  in  the  ionosphere  due  to  a  near  surface  underground  explo¬ 
sion.  In  this  case  the  explosion  was  taken  to  be  a  tamped  underground  nuclear  test  at  a  depth 
of  300  meters  with  a  seismic  body  wave  magnitude  near  5.  (Much  smaller  industrial  explo¬ 
sions,  very  near  or  at  the  earth’s  surface  ,  would  typically  produce  comparable  or  even  larger 
signals).  Figure  (27)  shows  electron  density  fluctuations  for  other  source  pressure  histories 
and  indicates  the  sensitivity  of  the  perturbation  form  to  source  character. 

4.  Modeling  High  Frequency  Seismic  Noise:  Atmospheric  Sources 

The  nature  of  high  frequency  seismic  noise  is  indicated  by  the  observed  noise  accelera¬ 
tion  power  shown  in  Figure  (28).  The  station  shown  (BAY)  is  near  the  former  Soviet  test  site 
in  Kazakh  and  is  typical  of  the  high  frequency  noise  seen  both  within  the  former  USSR  and 
elsewhere.  Three  components  of  ground  acceleration  on  the  surface  and  at  about  100  meters 
depth  are  shown.  Both  high  and  low  wind  level  seismic  noise  spectra  are  shown  on  each  plot. 
In  all  cases  the  high  frequency  seismic  noise  increases  with  high  wind  levels.  It  is  also 
apparent  that  the  acceleration  power  is  roughly  constant  over  the  band  from  about  1  Hz  to  30 
Hz.  Above  30  Hz  the  noise  acceleration  power  decreases  with  increasing  frequency,  particu¬ 
larly  in  the  bore-hole  at  100m.  depth. 

Given  the  rather  strong  dependence  of  the  noise  level  on  wind  velocity,  it  is  natural  to 
infer  that  atmospheric  coupling  at  the  earths  surface  is  an  important  means  of  excitation  of 
high  frequency  seismic  noise.  A  more  detailed  understanding  of  the  atmospheric  excitation  of 
seismic  noise  is  clearly  important  since  the  reduction  or  cancellation  of  this  noise  is  dependent 
on  an  understanding  of  its  origins,  mode  of  excitation  and  propagation  within  the  medium. 
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Figure  28.  Observed  high  frequency  seismic  noise  acceleration  power  spectra  at  a  site  near 
the  test  site  in  Kazakh.  Two  wind  level  conditions  are  shown  for  sensors  at 
the  surface  and  at  a  depth  of  about  100m.  (Data  from  Berger  et.  al„  1989). 


In  order  to  investigate  the  proati*.  ion  of  seismic  noise  by  atmospheric  processes,  the 
atmospheric  modeling  programs  were  linked  with  the  linear  elastic  seismic  modeling  pro¬ 
grams.  The  lower  atmosphere,  composed  of  a  day-time  turbulence  boundary  layer  with  a 
heigh  of  2  km,  is  simulated  with  a  random  surface  topography.  Winds,  blowing  on  the  topog¬ 
raphy,  induce  upward  and  downward  flow  velocities.  Random  temperature  changes  in  space 
on  the  ground  surface  also  produce  flows  that  self-organize  into  plumes  that  coalesce  above 
the  boundary  layer  into  larger  scale  thermals.  Together  with  random  turbulence  in  the  boun¬ 
dary  layer,  these  flows  induce  pressure  and  velocity  fluctuations  along  the  ground  surface. 
These  effects  are  the  input  into  the  seismic  modeling  code  which  integrates  in  time  from  the 
top-most  surface  boundary. 

Figure  (29)  shows  results  of  modeling  atmospheric  variations  in  pressure  and  velocity 
due  to  turbulence  (top  row)  and  induced  seismic  effects  at  depth  produced  by  these  atmos¬ 
pheric  effects  at  a  particular  time  (bottom  row).  An  interesting  feature  of  the  seismic  noise, 
produced  by  the  essentially  random  fluctuations  in  the  atmosphere,  is  its  relative  organization 
into  spatial  zones  of  coherent  compression  or  dilatation  (for  example)  at  a  particular  time. 
This  self-organization  of  the  seismic  noise  field  into  coherent  spatial  zones  in  response  to  ran¬ 
dom  atmospheric  excitation  is  characteristic  of  the  results  of  this  modeling  and  suggests 
approachs  to  high  frequency  noise  minimization  using  arrays  that  take  advantage  of  these  pat¬ 
terns. 

Preliminary  results  indicate  that  the  seismic  noise  that  is  produced  decreases  in  amplitude 
with  depth  and,  as  shown  in  Figure  (30),  produces  a  seismic  velocity  spectrum  that  has  a 
trend  that  decreases  as  1/f  with  increasing  frequency,  in  the  range  from  about  1  to  50Hz. 
Below  about  40  meters  the  seismic  noise  appears  to  interact  in  such  a  manner  that  much 


59 


SEISMIC  ATMOSPHERIC 
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Figure  29.  CROSS-SECTIONS,  50km  by  50km,  of  atmospheric  and  set  mic  fields  due  to 
atmospheric  turbulence  generated  by  random  temperature  perturbation.,  decaying  inversely 
with  height  and  possessing  a  uniform  random  distribution.  Blue  colors  denote  negative  values 
with  red  denoting  positive  values.  The  four  cross-sections  on  the  left  are  height  vertical  and 
width  horizontal;  the  two  cross-sections  on  the  right  are  width  and  breadth. 
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smoother  variations  in  spatial  distributions  are  obtained  than  at  the  surface  and  with  associated 
decreasing  fluctuations  in  time.  Both  topography  and  winds  are  found  to  be  of  major  impor¬ 
tance  in  terms  of  amplitude  and  character  of  the  noise.  Prom  preliminary  results  it  can  be 
expected  that  time  of  day  will  also  be  important  due  to  the  change  of  the  turbulent  boundary 
layer  with  the  heating  of  the  Sun  and  its  temporal  dependence. 


Figure  30 .  Theoretically  predicted  seismic  noise  due  to  models  of  atmospheric  turbulence  at  the 
earth’s  free  surface. 


Comparison  of  the  modeling  results  in  Figure  (30)  with  the  observations  in  Figure  (28) 
shows  that  the  velocity  predictions  are  in  agreement  with  the  observations  in  that  both,  in  the 
mean,  show  particle  velocity  decrease  as  1/f  with  increasing  frequency  above  about  5  Hz. 
This  strongly  implies  that  atmospheric  sources  (turbulence)  is  the  major  source  of  high  fre¬ 
quency  seismic  noise. 
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5.  Summary  and  Conclusions 


Since  the  atmosphere  is  just  another  layer  of  the  planet  and  wave  phenomena  in  this 
region  is  coupled  to  that  in  the  solid  layers,  and  is  of  a  form  that  is  only  modified  from  that 
of  seismic  waves  by  the  differences  in  material  properties  and  a  different  approximation  of  the 
same  equations  of  motion,  it  seems  evident  that  the  boundary  between  the  atmosphere  and  the 
lithosphere  should  not  constitute  a  barrier  to  investigation  or  utilization  of  signal  information. 
To  ignore  the  atmosphere  is  to  neglect  an  important  source  of  data  bearing  on  present  moni¬ 
toring  problems. 

In  addition,  it  might  be  pointed  out  that  sensing  of  atmospheric/ionospheric  disturbances 
by  electromagnetic  sounding  methods  is  only  utilization  of  a  particular  sensing  system  and  can 
be  considered  to  be  like  the  use  of  a  seismic  transducer  to  record  near  surface  ground  motion. 
Therefore,  differences  in  the  proposed  study  of  atmospheric  waves  for  monitoring  and  seismic 
monitoring  research  are  not  as  great  as  might,  at  first,  be  thought. 

The  potential  benefits  of  joint  seismic/atmospheric  monitoring  are  numerous.  Direct 
recording  of  acoustic  or  gravity  waves  at  surface  stations  provides  another  means  of  locating  a 
source.  Of  most  importance,  however,  is  that  merely  recording  such  signals  above  noise  lev¬ 
els  indicates  a  shallow  explosive  source  and  not  a  decoupled  nuclear  test.  Likewise  the 
indirect  EM  sensing  of  the  much  larger  ionospheric  wave  disturbances  due  to  acoustic  and/or 
gravity  waves  should  be  diagnostic  of  source  depth  simply  by  the  inferred  size  of  the  signal, 
where  surface  or  very  near  surface  industrial  explosions  (of  any  type)  will  produce  much 
larger  signals  than  an  earthquake  or  a  tamped  or  decoupled  nuclear  test  of  comparable  size. 
Naturally,  deeper  earthquakes,  that  can  sometimes  be  confused  with  nuclear  tests  because  of 
lack  of  good  seismic  estimates  of  depth,  will  produce  no  atmospheric  effects  while  a  shallow 
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nuclear  test  should  produce  a  moderate  to  small  (but  detectable)  signal. 

Therefore,  the  EM  detected  ionospheric  waves  could  be  very  useful  in  rapidly  identifying 
industrial  explosions  and  isolating  (by  elimination)  potential  decoupled  or  tamped  nuclear  tests 
out  of  a  group  of  seismically  determined  shallow  events  with  explosion-like  characteristics. 
Further,  unrecognized  deeper  earthquakes  could  probably  be  eliminated  as  potential  explosions 
because  of  the  lack  of  a  detectable  ionospheric  disturbance,  although  careful  modeling  of 
seismic  source  types  at  various  depth  locations  will  be  required  before  this  could  be  applied 
with  confidence. 

It  therefore  appears  that  by  placing  radio  frequency  receivers  and  transmitters  in  a  distri¬ 
buted  network,  comparable  to  an  in-country  seismic  monitoring  network  of  about  30  stations, 
should  make  it  possible  to  provide  complete  atmospheric  monitoring  capability  on  a  continen¬ 
tal  scale.  In  principal  it  should  be  possible  to  identify  large  industrial  explosions,  with  high 
probability,  due  to  the  strong  EM  signal  shifts  observed  and  thereby  allow  combined  seismic 
and  atmospheric  monitoring  to  identify  low  yield  coupled  and  decoupled  nuclear  tests  in  a 
background  of  industrial  explosions  and  earthquakes. 

We  have  tested  and  modified  atmospheric  modeling  capabilities  to  include  the  most 
important  non-linear  effects,  in  particular  the  effects  of  sub-grid  scale  turbulence.  We  have 
made  good  progress  on  the  turbulence  phenomena  by  introducing  a  randomly  fluctuating  com¬ 
ponent  to  the  field  variables  (i.e„  pressure,  velocities,  temperature  and  density)  that  simulates 
sub-grid  level  turbulent  effects.  Results  are  encouraging  and  in  particular  give  stable  dynami¬ 
cal  solutions  to  test  problems  that  are  comparable  to  observations. 

The  seismic  investigations  have  focused  on  generation  and  testing  of  2  and  3D  finite 
difference  programs  that  incorporate  surface  topography,  medium  randomness  and  lateral  vari- 
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ability.  Adaptation  of  FFT  methods  coupled  with  moving  grids  have  been  successful  when 
tested  against  analytical  and  conventional  finite  difference  methods  and  can,  in  principle,  pro¬ 
vide  capabilities  for  predictions  of  wave  fields  in  heterogeneous  media  at  regional  distance 
ranges  using  moderate  sized  computers  (e.g.,  high  level  work-stations  such  as  the  Stellar 
Computer)  with  only  modest  core  size  requirements  (i.e.,  a  few  hundred  megabytes). 
Transmitting  grid  boundaries  have  also  been  developed  and  tested  with  success.  Analytical 
(modal)  theory  methods  are  being  developed  for  3D  laterally  varying  media,  to  be  used  along 
with  the  2D  theory  already  developed. 

Results  of  modeling  studies  and  their  comparisons  with  observed  data  have  shown  that: 

(1)  The  atmospheric-ionospheric  modeling  predictions  of  electron  density  fluctuations  from 
the  non-linear  "gravity  wave"  have  large  amplitudes  and  wave  forms  of  distinctive  form 
suggesting  easy  detection  by  EM  sounding  methods. 

(2)  Coupled  atmospheric-seismic  modeling  indicates  that  atmospheric  turbulence,  simulated 
by  random  fluctuations  in  state  variables  near  the  free  surface,  produces  high  frequency 
seismic  noise  with  spectral  character  close  to  that  observed;  that  is  with  a  velocity  ampli¬ 
tude  spectrum  varying  as  1/f  as  a  function  of  frequency  above  1Hz. 

(3)  Seismic  wave  field  modeling  in  complex  structural  models,  incorporating  rough,  near 
source  topography  and  fine  scale  randomized  layering,  produces  seismic  synthetics  hav¬ 
ing  the  complex  character  of  observed  seismograms.  Simple  vertical  and  horizontal  ran¬ 
domization  can  only  be  an  adequate  representation  of  earth  structure,  and  the  seismic 
wave  fields  observed,  within  a  few  tens  of  wavelengths  from  the  source.  Beyond  that 
distance  the  effects  of  strong  shallow  lateral  variations  in  both  average  and  random 
characteristics  of  the  earth’s  structure  produce  effects  in  the  observed  wave  field  that 
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become  of  first  order  and  therefore  important,  so  that  accounting  for  large  scale  lateral 
variations  is  necessary  to  explain  observed  seismic  wave  fields  in  the  regional  distance 
range. 

(4)  Preliminary  studies  of  particular  complex  seismic  wave  types,  such  as  Lg  and  Rg,  indicate 
that  only  .arge  scale  lateral  variations  in  structure,  in  combination  with  both  vertical  and 
lateral  randomization  can  explain  the  wave  forms  and  attenuation  characteristics 
observed. 
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